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ABSTRACT : We use a modified Thomas-Fermi approximation to estimate analytically the 
critical velocity for the formation of vortices in harmonically trapped BEG. We compare this ana- 
lytical estimate to numerical calculations and to recent experiments on trapped alkali condensates. 



I. Introduction : The experimental realization of Bose-Einstein condensation (BEC) in cold 
trapped alkali atoms has stimulated great interest in weakly interacting inhomogeneous quantum 
gases. The alkali condensates are excellent systems for the study of quantum fluids, because (i) they 
are dilute, and hence admit an accurate description by mean field theory; (ii) system parameters 
such as density, temperature, and trapping potential are under fine experimental control; and (iii) 
sensitive diagnostics (mostly optical) have been developed to probe condensate behavior. In partic- 
ular, alkali BEC may serve as a powerful laboratory to study superfluidity - both long-established 
questions about its breakdown, and newer questions about the effects of flnite size and spatial 
inhomogeneity. 

It is well-known theoretically that collective excitations in a weakly interacting BEC determine 
a superfluid critical velocity, below which the condensate flow relative to a perturbing object or 
potential is without drag. This critical velocity is given by the Landau criterion [1] : 



where E[I) is the energy of a condensate collective excitation of linear momentum (or impulse) 
/. For liquid ^He, the breakdown of superfluidity is observed experimentally to occur at a critical 
velocity that is much smaller than that due to the excitation of phonons or rotons [2]. Beginning 
with Feynman more than 40 years ago [3], it has long been thought that the critical velocity for 
superfluid "^He is set by the creation of complex vortex structures (e.g., pairs, rings, and loops) 
depending on system geometry, temperature, etc. However, it has not been possible to fully verify 
this basic hypothesis because of the strong interactions in a liquid, which complicate quantitative 
comparison of superfluid ^He experiments [4] with theory [5], as well as the limited ability to vary 
the liquid density. Alkali gas BEC may provide an experimental system to test quantitatively the 
link between the superfluid critical velocity and the creation of vortices. 

In a landmark recent experiment [6], Ketterle and co-workers observed a critical velocity for 
the excitation of an alkali BEC when a perturbing potential (a blue-detuned laser beam) was moved 
through the trapped quantum fluid. These measurements are in qualitative agreement with both the 
well-known analytical calculation for a homogeneous weakly-interacting BEC in a channel of flnite 
diameter [2,3,7] (providing a Vc that is a factor of ~ 2 lower than experiment), and recent numerical 
calculations [8-11] based on the Gross-Pitaevskii equation [12] (providing a Vc that is a factor of ~ 
2 higher than experiment). 

Detailed understanding of the breakdown of superfluidity in Bose condensates will likely require 
further experimentation (e.g., to vary the trapped alkali BEC geometry and density, and to create 
greater spatial symmetry in the BEC perturbation), as well as additional theoretical work - both 
numerical and analytical. To this end, we report in this paper an analytical calculation of the critical 
velocity for vortex pair production in a harmonically-trapped dilute Bose condensate. We employ a 
modiflcation of the usual Thomas-Fermi (TF) approximation to the Bogoliubov mean-fleld theory 
[13], treating the BEC like a fluid with an exotic equation of state. We include leading kinetic 
energy terms caused by the two counter-rotating vortices, and use these kinetic energy terms as an 
effective potential for the BEC. We then compute, approximately, the energy E{I) and impulse / of 
the vortex pair, and determine the critical velocity from the Landau criterion. 

As a step toward checking this modifled TF limit for BEC vortices, we also compute the Bo- 
goliubov many-body "wavefunction" ('0) for a condensate with a single central vortex, and then 
calculate E, I, and Vc for this case. Since our calculations are based on an effective hydrodynamic 
model of BEC, the vortices that we consider do not, strictly speaking, have quantized angular mo- 
mentum in units of h. However, we find that the modified TF prediction agrees quite well with a 
vortex solution computed numerically. 
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We believe the advantages of the analytical method presented here are (i) its simplicity; and 
(ii) that it works well for the limit of large inter-particle interaction, which is the limit largely 
being explored by current experiments. In addition, the modified Thomas-Fermi analysis suggests 
the manner in which fraction of the speed of sound (cb), decreases as the inter-particle 

interaction increases. 

In section II below we introduce both the model used in our calculation, and the modified 
Thomas-Fermi approximation. In section III we apply this model and approximation to the creation 
of vortex pairs near the center of a harmonically-trapped BEC, and then compare the analytically 
estimated critical velocity with the recent measurements by the MIT group [6]. In section IV we 
show that a modified TF analytical calculation of Vc for a single central vortex in a trapped BEC 
agrees well with numerical simulation. Finally, in section V we summarize the implications of this 
approach. 

II. Model of Trapped BEC and the Modified Thomas-Fermi Limit : We restrict our 
discussion to the dynamics of a harmonically-trapped single-component Bose condensate at T = 0, 
which is well described by the Gross-Pitaevskii (GP) equation [12]: 
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where the field i/^', the diagonal part of the multiparticle wave function, serves as the condensate 
order parameter, and is referred to as "the wavefunction of the condensate"; M is the atomic mass; 
li' is the BEC chemical potential; and U\'ip'\'^ is an effective nonlinear potential arising from the 
atomic interactions (due to s-wave binary atomic collisions). Here U = , where is the 

atomic s-wave scattering length. For simplicity, we assume the trap has cylindrical symmetry and 
take R"^ = x'^ + y'^ — \f l"^ , with force constant K and no trapping potential along the symmetry (z) 
axis. Note that we solve the above GP equation subject to the constraint that the total number of 
atoms N is fixed, 

N = J d^rdz . (2.2) 

For ease of calculation we use scaled harmonic oscillator units (h.o.u.) in which the units of length, 

time, and energy are sj^^^^ , ^, and hoo, respectively, where uj = -y/^ is the trap angular frequency. 
Also, we reduce the GP equation to two dimensions by assuming the wavefunction is separable: 
= ■?/;(r )^{z). Thus, adopting the notation of Jackson et al. [10], we find 
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where we have scaled ip so that 
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Solving Eq. (2.3) subject to Eq. (2.4) fixes fi, the dimensionless chemical potential. The dimension- 
less interaction parameter is 

C=?^ (2.5) 

where 



is the classical turning-point width of the condensate in the harmonic trap in the C — > limit (equal 

to the unit of length in h.o.u.), and Q is the trap aspect ratio given approximately by C = |<^(o)P • 
For large C - which corresponds to most current experiments - the solution to the GP equation 
is well approximated by the Thomas-Fermi (TF) limit [13]. In this limit one neglects the derivative 
term in Eq. (2.3), yielding 




m) = ]l ^ ^'^ and , (2.7) 

where the second relation comes from the normalization condition, Eq. (2.4). (See also Eq. (3.4) 
and discussion below.) Eqs. (2.7) will be used extensively in this paper as the no-vortex solution 

that we compare to the vortex solutions. For reference, note that in the TF limit (i.e., the long 
wavelength limit), the local speed of (Bogoliubov) sound in the condensate is given in h.o.u. by 
cb — a/2C|'0P (see for example Ref. [10]). 

The many-body wavefunction of a vortex in the condensate will have a spatially-dependent 
phase [see Eq. (3.1) below], and a spatially-dependent amplitude which we refer to as the 'envelope 
function'. Let (p{f ) and A{f ) be the phase and envelope functions (both real) parameterizing 
this wavefunction via '0 = Ae^'^. Using this parameterization in Eq. (2.3), and equating real and 
imaginary parts, we find: 

-^"^ A + A{i<pf + ^A + CA^ ^ fiA (2.8) 

- V' - 2v^ • V0 = . (2.9) 

One may think of Eq. (2.8) as an equation of hydrostatic equilibrium and Eq. (2.9) as a conti- 
nuity equation. What we call the modified TF approximation consists of first solving the continuity 
equation [Eq. (2.9)] for (p{f) , and then using that solution in Eq. (2.8) to determine A{r), neglecting 
the Sj'^A term. Thus, in this approximation the envelope function A{f) is the solution of a simple 
algebraic equation. 

Once we have normalized the vortex wavefunction via Eq. (2.4), we compute the energy and 
integrated impulse. The energy functional is 

E = jd?r [yr ■ vV' + + '^{r^f] ■ (2.10) 

Using — Ae^'^ with A and 4> real, one finds, 

E = y dV [ivAf + {A^cpf + ^ + ^A^] ■ (2.11) 

The local momentum density is P = ^((vV'*) V' ~ V'*V'0); which is used to compute the total 
impulse: 

/ = J d^r \P\ 



= hj dVyl^lv*/*! . (2.12) 



We can then estimate the critical velocity for vortex pair creation by applying the Landau criterion 
[Eq. (1.1)], 

^vortex -E'no vortex /<-, -, on 

Vc = — — . (2.13) 
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In the next section we describe the details of an analytical computation of this ratio for a vortex 
pair near the center of the trap. 



III. Analytical Calculation of the Vortex Pair Critical Velocity : We consider a dilute BEC 
containing a pair of counter-rotating vortices with opposite charge or vorticity ±n. The vortices 
arc assumed to be located symmetrically about the trap center, with the cores a distance d apart. 
We apply the modified Thomas-Fermi (TF) method to the hydrostatic equilibrium expression [Eq. 
(2.8)] by assuming that is negligible. In addition, we adopt the ansatz that the two vortices are 
far enough from each other, but near enough to the trap center, that the total phase advance about 
the vortex pair is the simple sum of the phase advances of the two vortices viewed individually. This 
ansatz is akin to the analytical approximation developed by Feynman [3] for homogeneous BEC, 
and also used by Fetter in a similar context [14]. Practically, the ansatz is equivalent to requiring 
< d^ < fj, in dimensionless harmonic oscillator units (h.o.u.). Operationally it means that the 

wavefunction's phase (/>(r) satisfies = everywhere outside the vortex cores, which implies from 
the continuity expression [Eq. (2.9)] that ' V<^ must be zero in this region. Our solution for 
A{r) in this case satisfies the continuity equation near the vortex cores, and also far away from the 
vortices, since radial gradients of the phase vanish as In addition, the spatial integral of the 
continuity equation vanishes identically, and deviations due to the approximations described above 
are asymptotically bounded. 

We have found that neglect of the sy^A term in Eq. (2.8) leads to the main systematic errors in 
the analytical calculation described here. However, the deviations are small for parameters typical 
of current experiments. 

With the approach outlined above, the phase of the vortex pair wavefunction is 
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(3.1) 



where R and 9 are the usual polar coordinates as measured from the trap center. Employing the 
modified TF method, we find the condensate wavefunction envelope to be 

2 _R^ « \ 

~ Cy" 4 i?2d2 cos2 9+{R^- d^ /4)2 J ^ 

where the last term is the leading kinetic energy contribution, due to the |v0P term in Eq. (2.8). 
Note that this kinetic energy contribution is simply j^p^qyi, where ri and r2 are the respective 
distances measured from the vortex cores to the point {R,9). 

The above calculation scheme for A{r) breaks down very close to the vortex cores (small ri or 
r2) and also at large R, since the right hand side of Eq. (2.2) becomes negative. Considering these 
regions to be excluded complicates the analytic evaluation of energy and impulse [Eqs. (2.11) and 
(2.12)] precisely where the TF approximation fails. To circumvent this calculational difficulty we 
adopt a regulated expression for A{r): 

CX" 4 i?2d2cos2^+(i?2_rf2/4)2 + ^^ ^'^•'^^ 

^2 2 

where e = This regulated expression is justified by the observation that for vortex pairs not 

too far from the trap center (d^ < //), the contribution to A{f) from the kinetic energy term |V0P 
is never larger than ji. The regulation markedly affects the wavefunction envelope near the vortex 



cores, but has a small effect on the calculated energy, impulse, and critical velocities. (This is shown 
for a trap-centered, single vortex in section IV below.) In sum, the use of the regulated expression 
for A{f) is an important practical step in our analytical calculation, because it allows us to perform 
the spatial integrals for Ey^rtex and lyortex over the entire plane, rather than excluding regions near 
the vortex cores. 

We begin the integrations by normalizing the condensate wavefunction for the cases of no 
vortices and a single vortex pair. Referring to Eq. (2.7), the normalization for the no vortex case 
yields 

/ (J 

= /^o = Y 2^ , (3.4) 

where we have performed the integration in Eq. (2.4) out to a maximum radius, Rtf — 2^/1, 
i.e., the Thomas-Fermi condensate edge (the radius at which the condensate chemical potential is 
dominated by the trap potential in the Thomas-Fermi limit). Similarly, inserting Eq. (3.3) in Eq. 
(2.4), and noting that \ip\'^ = A^, we determine the normalization condition for the vortex pair 
wavefunction: 
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(3.5) 



where we have again performed the spatial integration out to the Thomas-Fermi condensate edge. 



which for large C is Rtf = 2^//U — ^^^/2 + • • •• Note that Eq. (3.5) indicates that ii increases with 
the vortex charge, n, as expected. (See Ref. [15] for details of the calculation of these normalization 
conditions.) 

Next, we compute the condensate energy. For the case of no vortex, we insert Eq. (2.7) in Eq. 
(2.10), and find 

1# ■ (-) 



For a vortex pair of charge ±n, we insert in Eq. (2.11) both the vortex pair's phase given in Eq. 
(3.1) and the regulated expression for A{f) given in Eq. (3.3). Performing the spatial integration 
again out to Rtf (see [15] for details), we find the energy for a vortex pair to be: 
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(3.7) 



Using Eqs. (3.6) and (3.7), we compute the energy difference between a trapped condensate with a 
symmetric pair of unit-charged, counter-rotating (n = 1) vortices, and the no vortex ground state 
(n = 0). To leading order in the limit of large C, we find 



E^-Eo = ^Hd'^io + 1) + . . . 



(3.8) 



Note that the energy difference vanishes in the d — > limit, as expected. Also, the energy difference 
varies as the natural logarithm of the distance between the vortices (at large separation), which 
resembles results from earlier studies of vortices in homogeneous BEC [2]. 

Next, we calculate the impulse In for the condensate with a symmetric pair of counter-rotating 
vortices with charge ±n. (With no vortices, the trapped condensate impulse is zero.) Using Eq. 
(3.1) we find 

IV0I = — r • (3.9) 

' ' [i?2(i2cos2 0+(i?2-d2/4)2]i 



We insert this formula and the regulated expression for A2(r) [Eq. (3.3)] into Eq. (2.12) to determine 
the vortex pair impulse. We find the largest contribution to the spatial impulse integral comes at 
large R, and is thus dependent on the condensate size (set roughly by the Thomas-Fermi condensate 
edge, Rtf) as well as the vortex pair charge magnitude (n) and separation {d). Evaluating the 
integral out to Rtf, and assuming large C, we find the vortex pair impulse to be 

/„ = -^ln(-^) + ... (3.10) 

to leading order (sec [15] for details). Note that like the energy, the impulse for vortex pair creation 
vanishes in the small d limit, and is proportional to the vortex charge, as expected. Also, since we are 
using a macroscopic hydrodynamic model in our calculations, the quantization of angular momentum 
of the BEC wavefunction does not lead to a simple quantization condition on the integrated impulse. 

Assembling the results for the energy and impulse leads directly to an estimate of the critical 
velocity at which symmetrically- placed, oppositely-charged vortex pairs, at relative separation d, 
can form near the center of a Bose condensate trapped harmonically in two dimensions. In the large 
C limit by using Eq. (2.13) we find 

21n(^ (3.11) 



where again, /xq = y 5^- These expressions indicate a weak dependence (via jiQ and C within the 

logarithms) of Vc on the BEC non-linear self-energy. Recall that the derivation of Eq. (3.11) assumes 
that 1/ n < d^ < n m. dimensionless harmonic oscillator units (h.o.u.). The first inequality is the 
requirement that the vortices be at least several healing lengths away from each other, and the 
second inequality is the requirement that they are not too far from the trap center. 

These conditions seem to be largely met in the recent MIT critical velocity experiment using a 
trapped sodium BEC [6] , although this experiment and our model have different spatial symmetries 
(i.e., geometries). Onset of energy dissipation in the MIT condensate was observed at a critical 
velocity when moving a repulsive barrier (a blue-detuned laser beam) through the middle of the 
cold atom cloud. The laser beam was directed along the smaller, radial direction of the cigar-shaped 
condensate, which had Thomas-Fermi diameters of 45 /xm and 150 /xm in the radial and axial 
directions, respectively. The laser barrier was moved back-and-forth along the condensate's axial 
direction (perpendicular to the axis of the laser beam) at a constant speed, and energy dissipation was 
determined from measured changes in the condensate fraction. Thus the experiment was explicitly 
a three-dimensional system with anisotropic trapping in the plane perpendicular to the laser beam 
axis. 

Ignoring the different geometry of our two-dimensional, isotropically-trapped BEC model, we 
use the analytical calculation outlined above to estimate the critical velocity for vortex pair creation 
in the MIT condensate. We assume the vortex cores are parallel to the axis of the blue-detuned laser 
beam, and are separated by a distance d equal to the laser beam diameter (13 //m). In the central 
region of the MIT condensate, the number of atoms per unit length parallel to the laser beam axis 
{til) was about 3 x 10^ cm"-*^. Using 2.9 nm for the sodium s-wave scattering length (a^), we find 
from Eq. (2.5) that C ~ SnasTiL ~ 25,000 (in h.o.u.) for the MIT BEC experiment. 

Using these values in Eq. (3.11), we estimate that the onset of vortex pair formation, and hence 
energy dissipation, occurs at a critical velocity Vc — 0.4 mm/s, a factor of about four below the 
measured value [6]. (Note: this estimated Vc varies by ~ 30% over the range of trap frequencies 
present in the MIT experiment - from 18 Hz in the axial direction to 65 Hz in the radial direction.) 
We defer further interpretation of this estimate of Vc to the conclusion. 



IV. A Test of the Modified Thomas-Fermi Limit : As a test of the analytical calculation 
of vortex critical velocity using the modified Thomas-Fermi (TF) approximation, we consider a 
simple, symmetric example: a single vortex of charge n at the center of an isotropic BEC that 
is harmonically-trapped in two dimensions. In this case, the cylindrical symmetry of the system 
requires the phase of the condensate to advance proportional to the polar angle 0, and hence Eq. 
(2.9) limits the envelope function A{r) to be a function oi R = \f\ alone. Also, the single- valuedness 
of the condensate wavefunction constrains the phase to be = nO, where n is an integer (the vortex 
charge). Therefore, Eq. (2.8) reduces to a second order ordinary differential equation for A{R) that 
can be integrated numerically, without making the modified Thomas-Fermi approximation, which 
can then be used to test the modified TF analytical calculation. We employed an adaptive mesh 
relaxation algorithm to perform the numerical integration of Eq. (2.8) for both n = and n — 1, 
using the value C — 200,000 for the coefficient of the non-linear term throughout. The computed 
A{R) for n = is shown in Figure 1 and compared to the corresponding A{R) calculated analytically 
using the traditional TF approximation. 
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Fig. 1: Calculated envelope function A{R) for a condensate 

with no vortex (n = Q). Numerical and TF analytical 
calculations give identical results on the scale of this graph. 



Figure 2 provides a comparison of the numerical and modified TF analytical solutions for A[R) 
with a single central vortex of charge n = 1. 
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Fig. 2: Comparison of numerical and modified TF envelope functions A{R) 
for a condensate with a single central vortex of unit charge {n = 1). 
(a) Full radial extent of A{R): differences not perceptible on this scale, 
(h) A{R) near the trap center. Dashed line is modified TF calculation. 

In the modified TF limit, the analytical envelope function for a single central vortex of charge 

n is 

MR) = ^ ^ , m = ne . (4.1) 

This expression for A{R) shows that including the leading kinetic energy term due to the vortex in 
the analytical calculation is equivalent to an effective potential barrier near the vortex core caused 
by the angular momentum of the condensate. In this approximation the condensate extends over 
the annulus Rq < R < Rmaxi outside of which A vanishes. From Eq. (4.1) we see that 

Rl = 2ii- _ ^2 and i?^^^ = 2/1 + 2^p? - n? . (4.2) 

Hence in the modified TF limit, the vortex core radius (in h.o.u.) is effectively n/ at large ii. On 
general grounds, we expect the radius of the vortex core to be of the order of the condensate healing 
length, ^ = g^^^ , where is the s-wave scattering length and no is the local density. Thus the 

modified Thomas-Fermi limit reproduces the expected scaling of vortex core size, since ^ ~ in 
dimensionless units. Note also that at large R (ignoring the trap potential) the effect of the vortex 
on the condensate envelope falls off as 1/R^ in the modified TF approximation, which is identical to 
the asymptotic behavior of the trial vortex wavefunction in homogeneous BEG used by Fetter [14]. 

The modified TF analytical calculation of the energy of a single central vortex of charge n is 
found by using Eq. (4.1) in Eq. (2.11) (ignoring the term) and integrating over the annulus 
given in Eq. (4.2): 



^ 47r(u^ — n^)2 

En = ^^3^ ' . (4.3) 

[Compare with Eq. (3.6).] Although not obvious from this equation. En increases as a function of n 
due to the dependence of /x on n through the normalization condition J d^r \ip\'^ — 1. In particular, 
we find this normalization requires 



^ = ,V^?^-^4 ^'^^gIg ) . (4.4) 
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hence causing Rmax (Eq. (4.2)) to increase with n. Next, using Eq. (4.1) in the analytical calculation 
of the impulse [Eq. (2.12)] of a single central vortex, we find 

^ = ^(/^-^)^ • ^^-^^ 

For the no-vortex case {n = 0) and C = 200, 000, numerical integration of Eq. (2.3) yields the total 
condensate energy: Eq = 118.965. [All numbers here are in harmonic oscillator units (h.o.u.).] With 

the analytical method we find Eq = = 118.94, where ijq = \J~%^ = 178.4. For a single central 
vortex with n — 1 and C — 200,000, the numerical solutions give Ei — 118.991 and Ii — 0.0996, 
whereas the modified TF analytical calculation yields Ei = 118.97 and Ii = 0.0989. Using these 
values, we find Vc = = 0.26 (numerical) and 0.30 (analytical). [Recall that the impulse integral 
is identically zero for the no vortex solution.] This reasonable agreement suggests that the modified 
Thomas-Fermi approximation in the large C limit should enable analytical calculations of 10-20% 
accuracy for vortex critical velocities in inhomogeneous BEC, provided the calculational model and 
experiment have the same spatial symmetries. 

IV: Conclusion : In this paper we report an analytical calculation of the critical velocity for 
creation of a pair of counter-rotating vortices in harmonically-trapped, dilute Bose condensates. 
Excitation of vortices may set the lowest limit for a superfluid critical velocity in BEC. For the 
analytical calculation presented here, we developed a modified Thomas-Fermi approximation to the 
standard mean-field theory [13]. This approximation is appropriate to the limit of large intcrparticle 
interactions, and thus is relevant to most current experiments. The approximation consists of two 
steps, described above in Section II: (i) solve a continuity equation for the phase of the condensate 
wavefunction; and then (ii) insert this solution for the phase in an equation of hydrostatic equilibrium 
and solve for the condensate envelope function, ignoring the second order spatial derivative of the 
envelope function. In essence, we include contributions to the condensate kinetic energy that come 
from gradients of the phase (a signature of vortices), but assume that vortex- induced gradients in the 
condensate envelope function are negligible. Alternatively, one can think of the modified Thomas- 
Fermi approximation as equivalent to treating BEC like a fluid with an exotic equation of state, 
with the leading vortex kinetic energy terms acting as an effective potential. Once the envelope 
function is determined, we compute approximately the energy and integrated linear momentum 
(or impulse) of the vortex pair, and determine the critical velocity from the Landau criterion, 
Vc = miniener gy / impulse) . 

As described in Section III, we find qualitative agreement between our analytical calculation and 
the BEC excitation critical velocity recently measured by Ketterle and co-workers at MIT [6] . This 
result is not surprising given the obvious difference in symmetry between our model and the MIT 
critical velocity experiment: we assume the vortex core axes to be parallel to the axis of cylindrical 
symmetry in a two-dimensional trap confining the condensate; whereas, the experimental excitation 
would likely create vortices with core axes perpendicular to the finite cylindrical axis in the actual, 
quasi-two-dimensional trap. In addition, vortices produced at velocities substantially below the 
observed critical velocity may collectively dissipate energy at a rate comparable to the decay rate 
of the unstirred condensate itself, making it difficult to identify Vc experimentally (i.e., the MIT 
experiment may only place an upper limit on Vc)- We expect that detailed understanding of the 
role of vortex creation in the breakdown of superfluidity in BEC will require further experimental 
development to enable more symmetric excitation of the condensate, as well as theoretical treatments 
(both numerical and analytical) of more realistic systems. Nevertheless, as shown in Section IV, the 
modified Thomas-Fermi analytical method is in reasonable agreement with numerical calculations 



based on the full GP equation of the critical velocity for a vortex with high symmetry: specifically, 
a single vortex in the center of a two-dimensional harmonic trap. 

We note two important issues that are not included in the analysis presented in this paper: 

(i) the breakdown of the hydrodynamic approximation at the edges of the trapped condensate; and 

(ii) the dynamics of vortex excitation in the presence of nominally smooth trapping and perturbing 
potentials. Further theoretical and experimental work is necessary on these problems. 

We conclude by emphasizing that the modified Thomas-Fermi approximation presented in this 
paper is straightforward to implement. In addition to enabling an analytical calculation of vortex 
critical velocity (for more details see Ref. [15]) this approximation should have other applications 
in the physics of trapped BEG. 

We thank E. Timmermans for helpful discussions. 
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